void flow_study(){

	TH1F *errevt2 = new TH1F("h1","Event plane 2 error",200,-2,2);
        TH1F *errevt3 = new TH1F("h2","Event plane 3 error",200,-2,2);
        TH1F *errv2 = new TH1F("h3","v2 error",200,-2,2);
        TH1F *errv22 = new TH1F("h4","v2 error2",200,-2,2);


	TRandom *rndm = new TRandom3(0);
        int nevent = 1000;
        int nparticle = 1500;
	TF1 phidist("phidist", " [0] + [1]*TMath::Cos(x-[2]) + [3]*TMath::Cos(2*(x-[4])) + [5]*TMath::Cos(3*(x-[6])) ", 0, 2*TMath::Pi());
        const double a0 = 10;



      for (int ievent=0;ievent<nevent;ievent++){
        vector<double> tracks;
	double An = 0;
        double Bn = 0;	
        double An_3 =0;
        double Bn_3 =0;
        double v2_measure2=0;
	phidist.SetParameter(0, a0); // const
	phidist.SetParameter(1, 0); // direct flow value
	phidist.SetParameter(2, 0); // ingore 1th order event plane
	phidist.SetParameter(3, 0.2); // v2
	double evtplane2 = rndm->Uniform(0, 2*TMath::Pi()/2); 
	phidist.SetParameter(4,evtplane2); // event plane (will be change later)
	phidist.SetParameter(5, 0.2); // v3
	double evtplane3 = rndm->Uniform(0, 2*TMath::Pi()/3); 
	phidist.SetParameter(6,evtplane3 ); // evnet plane for 3th order
        double integral = phidist.Integral(0,2*TMath::Pi());
	for (int iparticle=0; iparticle<nparticle; iparticle++){
		tracks.push_back(phidist.GetRandom());  // filling particle
	}  // iparticle loop	
        for (int i=0;i<nparticle;i++) {

	An += TMath::Cos(2*tracks[i]);
        Bn += TMath::Sin(2*tracks[i]);
        An_3 += TMath::Cos(3*tracks[i]);
        Bn_3 += TMath::Sin(3*tracks[i]);

	} //i loop

//        cout << "An is : " << An << endl;
//        cout << "Bn is : " << Bn << endl;

        double evt2_measure = (TMath::ATan2(Bn,An))/2;
        double v2_measure = (TMath::Sqrt(An*An + Bn*Bn))/nparticle*2*a0;
        double evt3_measure = (TMath::ATan2(Bn_3,An_3))/3;

        for (int i=0;i<nparticle;i++) {
        v2_measure2 += TMath::Cos(2*(tracks[i] - evt2_measure))/nparticle*2*a0;
        } //i loop

//        cout << v2_measure << endl;

        errevt2->Fill(evt2_measure - evtplane2);
        errevt3->Fill(evt3_measure - evtplane3);
        errv2->Fill(v2_measure - 0.2);
        errv22->Fill(v2_measure2- 0.2);

	} // ievent loop
        TCanvas *c1 = new TCanvas();
        c1->Divide(2,2);
        c1->cd(1);
        errevt2->Draw();
        c1->cd(2);
        errv2->Draw();
	c1->cd(3);
        errevt3->Draw();
        c1->cd(4);
        errv22->Draw();


//	1) Evt - evtplane (2, 3th order)
//	2) v2 - 0.2
//      issues :  extracting v2 is something wrong!
//  
}
